Triangularity and Dipole Asymmetry in Heavy Ion Collisions 



Derek Teaney and Li Yan 

Department of Physics and Astronomy, 
Stony Brook University, Stony Brook, 
New York 11794-3800, United States 
(Dated: October 12, 2010) 

Abstract 

We introduce a cumulant expansion to parameterize possible initial conditions in relativistic 
heavy ion collisions. We show that the cumulant expansion converges and that it can system- 
atically reproduce the results of Glauber type initial conditions. At third order in the gradient 
expansion, the cumulants characterize the triangularity (r 3 cos 3(0 — 1^3,3)) and the dipole asymme- 
try (r 3 cos(0 — ^1,3)) of the initial entropy distribution. We show that for mid-peripheral collisions 
the orientation angle of the dipole asymmetry -01,3 has a 20% preference out of plane. This leads to 
a small net v\ out of plane. In peripheral and mid-central collisions the orientation angles ipi^ and 
■03 5 3 are strongly correlated, but this correlation disappears towards central collisions. We study 
the ideal hydrodynamic response to these cumulants and determine the associated v\/e\ and ^3/63 
for a massless ideal gas equation of state. The space time development of v\ and i>3 is clarified 
with figures. These figures show that v\ and V3 develop towards the edge of the nucleus, and 
consequently the final spectra are more sensitive to the viscous dynamics of freezeout. The hydro- 
dynamic calculations for V3 are provisionally compared to Alver and Roland fit of STAR inclusive 
\q two particle correlation functions. Finally, we propose to measure the v\ associated with the dipole 

asymmetry and the correlations between ipi^ and ^3,3 by measuring a two particle correlation with 
respect to the participant plane, (cos(0 Q — 30^ + 2^>pp)). The hydrodynamic prediction for this 
correlation function is several times larger than a correlation currently measured by the STAR 
collaboration, (cos(0 a + (pp — 2^pp)}. This experimental measurement would provide convincing 
evidence for the hydrodynamic and geometric interpretation of two particle correlations at RHIC. 
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I. INTRODUCTION 



In a recent and significant paper B. Alver and G. Roland (AR) [lj provided the most 
compelling explanation to date for the striking features measured in two particle correla- 
tions at the Relativistic Heavy Ion Collider [2H7] . These features (which are described with 
picturesque names such as the "ridge" and "shoulder") are said to arise from the collective 
response to fluctuating initial conditions. Specifically, if the initial conditions are parame- 
terized with a quadrapole and triangular moment, the two particle correlations reflect the 
response of the nuclear medium to these anisotropies. The work of AR was motivated in 
part by event by event simulations of heavy ion collisions with ideal hydrodynamics which 
showed that the flow from fluctuating initial conditions can describe the general features of 
the measured two particle correlations [8|. The general idea that the curious correlations 
are due to a third harmonic in the flow profile was previously suggested by Sorensen In 
addition, many of the features of the observed two particle correlations were found in the 
AMPT model p^0HT2] . though the geometric nature of these correlations was not understood 
before the work of Alver and Roland. 

The hydrodynamic interpretation of the measured two particle correlations is important 
for several reasons. First, before this conclusion there was a significant correlation between 
the measured particles which was not understood. This confusion casted doubt on the hy- 
drodynamic interpretation of RHIC results and clouded the important conclusion that the 
shear viscosity to entropy ratio of QCD is of order ~ h/Aix near the phase transition [T3] . 
However, since the unusual two particle correlations are actually a prediction of hydrody- 
namics, the observation of these unusual features in the data validates hydrodynamics as an 
appropriate effective theory for heavy ion events and marginalizes other models. Further, 
once the hydrodynamic interpretation is adopted the measured correlations can be used to 
constrain the properties of the medium, e.g. the shear viscosity and the Equation of State 
(EOS). In particular the effect of viscosity was calculated in Refs. [HI [T3] which will be 
discussed more completely below. 

Motivated by these results, the current work will characterize the fluctuating initial con- 
ditions with a cumulant expansion. Instead of running hydrodynamics event to event, the 
linear response to specified cumulants can be calculated with ideal or viscous hydrodynamics. 
Subsequently, these response functions can be combined with a Glauber model for the event- 
by-event cumulants (and their correlations), and the combined result can be fairly compared 
to data. At third order in the gradient expansion, the initial condition is parameterized by 
a radial dependence to the dipole moment, (r 3 cos0), and the triangularity, (r 3 cos 30) . In 
Section III we will calculate (with ideal hydrodynamics) how the medium responds to these 
moments and illustrate how this response develops in space and time. Subsequently in Sec- 
tion IV we will compute the corresponding particle spectra Ui(pr) and v 3 {px) and study the 



sensitivity to certain model parameters related to freezeout. In Section VA we will make a 
comparison to V3A/V2A as extracted by Alver and Roland in their analysis of two particle 
correlations. We will also make definite predictions for the dipole asymmetry vi(p T ), which, 
if confirmed, would firmly establish the geometric nature of the two particle correlations. 
The comparison to data is not final as the effects of resonance decays, viscosity, and higher 
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cumulants have not been included, Nevertheless, the preliminary comparison will firmly tie 
the formal cumulant expansion outlined in this paper to the measured correlations. Finally, 
we will compare our calculations to the recent results of Refs. [HHT6] in Section V B 



II. CUMULANT EXPANSION AND HYDRODYNAMICS AT RHIC 

A. The initial conditions for ideal hydrodynamics 

In this paper we will use 2+1 dimensional boost invariant ideal hydrodynamics to simulate 
RHIC events [TBI [T71 IT8] . Briefly, in ideal hydrodynamics the stress tensor satisfies the 
constituent relation and the conservation laws: 

r = (e + ?(e))«Y + P(e)f , V M T M1/ = , (2.1) 

where e is the energy density, u M is the flow velocity, and the pressure V is specified by the 
EOS, V = V(e). We will work in flat space but with coordinates 

T = y/t 2 - Z 2 , Tj s = - log 

With the assumption of boost invariance, the hydrodynamic fields are independent of rj s and 
vP = 0. Using the constraint u^ L = —1, the independent fields which must be determined 
by solving the conservation laws are 

e(r,x), u x (r,x) u y (r,x), (2.2) 

where x denotes two dimensional vectors in the transverse plane. We will specify the initial 
conditions for the subsequent evolution in what follows. At the initial time r Q it is reasonable 
to assume that flow fields are small, u x ~ u y ~ 0. This leaves the initial energy density 
which must be specified e(r Q , x, y). We will specify the initial entropy density s(r , x, y) with 
a cumulant expansion and infer the initial energy density from the equation of state. 

A typical initial condition might be fairly complicated involving several structures. How- 
ever, the effect of the shear viscosity is to damp the highest Fourier modes. Thus, after 
damping the shortest wavelengths, the initial entropy distribution is approximately described 
by a Gaussian with average squared radius (r 2 ) and elliptic eccentricity e 2 as has tradition- 
ally been used to characterize heavy ion events [T7j. The damping of the highest Fourier 
modes is nicely seen in Fig. 1 of a recent preprint [J3]. The next paragraphs formalize this 
description and categorize corrections. 

B. Cumulants 

The Fourier transform of the entropy density for a given initial condition is 

J d 2 x e ik - x p(x) = p(k) , (2.3) 
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where p(x) = t s(t , x) / S to t and 5 to t = / r d 2 x s(r G , x) is the total entropy per space time 
rapidity. Since the highest Fourier modes are damped, we will expand the initial distribution 



in k. Expanding both sides of Eq. (2.3) with respect to k 



(ik a )(ik b ) , , 

p(k) = 1 + ik a p 1>a + K £ V afr + • • • , (2.4) 

we see that p(k) generates moments of the entropy distribution 

Pl,a = (X a ) , p 2 ,ab = (x a X b ) , (2.5) 

where the average is appropriately defined 

(...)= / d 2 ^)-- - • (2-6) 



Although we could classify the initial conditions with these moments, a cumulant expansion 
seems more natural since the average Glauber distribution is roughly Gaussian and the 
cumulants are translationally invariant. We therefore define W(k) 

exp{W(k)) = J d 2 xe ikx p{x), (2.7) 

and expand both sides in a fourier series 

W(k) = 1 + ik a W ha + ^(tk a )(tk a )W 2M + ... . (2.8) 

From this expansion we see that W(k) is the generating function of cumulants of the un- 
derlying distribution p(k) 

W hl = (xi) , W 2:ab = (x a x b ) - (x a ) (x b ) . (2.9) 

From now on we will shift the origin so that (x a ) = 0, and the distribution is approximately 
Gaussian to quadratic order 

p(k)=e^(-^k a k b W 2 ,a^ . (2.10) 

Higher order corrections in this expansion will correct the distribution away from the Gaus- 
sian. The tensor W^ ab is a reducible tensor and should be decomposed into irreducible 
components, 

1„. , /„. 1, 



W 2 , ab = -W 2 , cc 5 ab + [W 2 ,ab - -W 2>cc 8 ah ) . (2.11) 



We orient the x, y axes to the participant plane [19] where W 2tXy = 0. Then the irreducible 
moments are 

W 2taa =(x 2 + y 2 ) , (2.12) 

W 2>xx - l -W cc 5 xx = \{x 2 ~ V 2 ) ■ (2.13) 
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Clearly the irreducible components of the cumulant expansion are related to the traditional 
parameters of heavy ion physics: 

(x 2 + y 2 ), and e 2 ^ ^ ~f ] . (2.14) 

To write down corrections to these results it is more convenient and illustrative to use 
cylindrical tensors rather than Cartesian tensors. Appendix [A] develops this expansion in 
detail and only certain features will be summarized here. Appendix A 1 expands W(k) in a 
Fourier series: 

oo oo 

W(k) = W (k) + 2 W:(k) cos(0 fc ) + 2 ^ W s n (k) sin(n0 fc ) , (2.15) 

n=l n=l 

where k and 0& are the norm and azimuthal angle of the momentum vector. The W£' s (k) 
are also expanded in k to characterize the distribution at largest wavelength: 

W (k) = ^W , 2 (tk) 2 + 0(k 4 ) , (2.16a) 

W1(k) =Wl l + 0{k 3 ) , (2.16b) 
Wl(k) =Wl A + 0(k 3 ) , (2.16c) 

Wm =^Wl 2 {ik? + 0(k A ) , (2.16d) 

W°{k) =Wl 2 + 0{k 4 ) . (2.16e) 
After Appendix [A] we find that to order k 2 

(r 2 ) , (2.17) 

(2.18) 
(2.19) 

<r 2 cos(20)) , (2.20) 

(2.21) 

Here we have used translational invariance and rotational invariance (as in the Cartesian 
case) to eliminate W° l5 W/ l5 and W^ 2 . To third order in the gradient expansion the dipole 
terms W{{k) and W((k) are non-zero 

WZ{k)=^W 1)3 {ik) 3 + 0{k 5 ), W h3 =^(r 3 cos 0) , (2.22a) 

W;(k)=^W lt3 (ik) 3 + 0(k s ), H/ li3 =^<r 3 sin0) . (2.22b) 

Similarly, at third order in the gradient expansion there are terms proportional to cos (30) 

W£(k) =±W 3 , 3 (ik) 3 + 0(k 5 ) , W£ 3 =^ <r 3 cos(30)> , (2.22c) 

Wl{k) =±W 3 , 3 (tk) 3 + 0(k 5 ) W£ t3 =i <r 3 sin(30)> . (2.22d) 
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Once the fourier coefficients W nm are specified, the entropy distribution in space can be 



found with a fourier transform; see Eqs. 2.34 and 2.35 and the surrounding text for further 
discussion. 



C. A strategy for event by event hydrodynamics 

If the cumulants beyond second order are in some sense small, then the change in the 
hydrodynamic spectra due to a specified set of higher cumulants is linearly proportional to 
the deformation 



dSN 

d(j) p 



E 

n,m,{s,c} 



1 dSN 



c,s 
n,m 



where 



1 dSN 



Wi 



P J n,m,{s,c\ 



yy n,m 



(2.23) 



(2.24) 



^■$P . n,m,{c,s} 

is the normalized response to a given cumulant. If the non-linear interactions between the 
elliptic flow and the higher cumulants can be ignored (i.e. the elliptic flow is sufficiently 
small), then the background Gaussian is approximately radially symmetric and the response 
of the sin terms are related to the response of the cosine terms through a rotation. In this 
case, we are free to rotate our coordinate system by an angle ^> nm 



n 



atan2(^ m ,W^J + 



7T 



n 



(2.25) 



so that the sin terms vanish. In this rotated frame (which we will notate as W) the cumulants 

are 

W s 







and the spectrum can be written 



W. 
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E 



d5N 



(2.26) 



(2.27) 



n,m,{c} 



Thus, the assumption of a rotationally invariant background reduces the number of coeffi- 
cients by a factor of two. 

In this paper we will assume that all deformations from spherical are small including 
the elliptic flow. Thus, we will neglect the non-linear couplings between the elliptic flow 
and the triangular flow and the elliptic flow and the dipolar flow. We have investigated the 
influence of the ellipticity on the triangular and dipolar flow and our preliminary findings 
show that the effect of the elliptic flow on the triangular flow is small. A similar finding 
was reported in the very recent preprint by the Duke group [20] . However, the effect of the 
elliptic flow on the dipolar flow is non- negligible when the dipole angle is oriented in plane. 
This complication will be reported on in future work [2"T] . 
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The angle ^2,2 specifies the orientation of the participant plane typp, and the second 
order cumulant W22 determines the ellipticity 

_ (r 2 cos 2(0-^)) _ 4Wfc _ 
e 2 = , W PP = ^ 2>2 . (2.28) 

The participant plane angle ^ pp is distinct from the reaction plane angle which we denote 
with 

The third order cumulant W 33 describes the triangularity as introduced by Alver and 
Roland. These authors suggested a definition of the triangularity and orientation angle e AR 
and ip 3 R with a quadratic radial weight 



ef R = - (r2C ° S W 2) *£*))) , *An = l atan2((r2 sin(30)) ; {r 2 cos(30))) + I (2 29) 

We will abandon this analytically frustrated definition, and define the triangularity 63 and 
the associated angle with an r 3 weight 

(r 3 cos 3(0 -y>3, 3 )) 

^3 i3 =^atan2«r 3 sin 30) , <r 3 cos 30)) + ^ . (2.31) 

The difference between the r 2 and r 3 weight is captured by the response of the system to 
the fifth order cumulants, W§ 5 oc [(r 5 cos 30) — 4 (r 2 ) (r 3 cos 30)] . Recent studies of the 
response of the system to €5 (or in the current context) suggests that the response to 
these fifth order cumulants will be small [14]. 

The third order cumulant W£ 3 describes a dipole asymmetry and also appears to the 
same order in the gradient expansion. By analogy we define 6\ and ipx^ 

(r 3 cos(0-^, 3 )) _ 8^3 
ei= p) -1(H)' (2 ' 32) 

i[) 1)3 =atan2((r 3 sin 0) , (r 3 cos 0))+7T. (2.33) 

Estimates for these parameters and their correlations will be given in the next section. 



D. The dipole asymmetry and triangularity 



To get a feeling for the dipole asymmetry and triangularity we first record the explicit 
coordinate space expressions for a distribution with only triangularity 



S{X,T) OC 



1 + 



(r 3 ) e 3 



24 



d_ 

dx 



d V d 



and a distribution with only a dipole asymmetry 



s(x, T) oc 



1 + 



(r 3 ) ei 



d_ 

dx 



dy J dx 



d\ d 



(2.34) 



dy J dx 



(2.35) 
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FIG. 1: A schematic of an event with (a) net triangularity and (b) net dipole asymmetry. The 
triangularity produces a net v%(pt) and the dipole asymmetry produces a net v\(pr)- The cross in 
(b) indicates the center of entropy (analogous to the center of mass) and the large arrow indicates 
the orientation of the dipole. 

Here the orientation angles ^3 and ipi^ are set to zero. At large enough radius the deriva- 
tive terms become large and overwhelm the leading term making the distribution negative. 
This is an unavoidable consequence of truncating a cumulant expansion at any finite order. 
As explained in Appendix [X] we regulate these terms and adjust the overall constant to 
reproduce the total entropy in a central RHIC collision. Fig. [T^i and Fig. [T]d illustrate initial 
conditions with net triangularity and net dipole asymmetry respectively. The distribution 
with net triangularity leads to a v${pt) while the dipole asymmetry leads to a vi(pt)- 

To estimate these parameters and their correlations we have used the PHOBOS monte 
carlo Glauber code [22]. Fig. [2] shows the distribution of e%, 62 and e 3 as a function of the 
number of participants. We see that the dipole asymmetry is about a factor of two smaller 
than the triangularity but is not negligibly small. 

Fig. [3] shows the distribution of ipi^ and ^3 with respect to reaction plane at various 
impact parameters. We see that although ^3 3 is uncorrelated with respect to the reaction 
plane, ^1,3 shows an anti-correlation with respect to the reaction plane, which eventually 
disappears toward central collisions. 

More importantly, the angles ipx^ and ip^ are strongly correlated in mid central collisions 
(a similar observation was made recently by Staig and Shuryak [23] )• Fig- [4] shows the 
conditional probability distribution, i. e. 



The strong correlation may be explained physically as follows. When the dipole asymmetry 
is in plane then the triangular axis is at 7r/3, i.e. the point of the triangle is aligned with 
the dipole axis as exhibited in Fig. [5]^a). However, when the dipole axis is out of plane then 
the triangular axis is also out of plane as exhibited in Fig. [5tb). 




These correlations are a reflection of the almond shape geometry and their general form 
can be established by symmetry arguments. First, since the probability of finding a dipole 
asymmetry in a given quadrant of the ellipse is the same for every quadrant, the probability 



-P(^3,3 |^i,3, *k) = The probability of ^3,3 given ^1,3 and ^ R . 
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FIG. 2: Size of the moments e±, €2 and €3 as a function of the number of participants. The points 
indicate the average value of ((e n )) and the errorbars indicate the variance of e n at fixed -/V part . 
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FIG. 3: Distribution of the angles -01,3 and ^3,3 with respect to the reaction plane for three 
different impact parameters. 



distributions dP/d(?/>i,3 — ^r) must only involve even cosine terms 

dP 1 



d^ 



1,3 



2tt 



(1-2Acos2(Vi, 3 -*«)) + ■■■) • 



(2.36) 



The sign has been chosen so that a positive A coefficient describes the out-of plane preference 
seen in Fig. |3j The coefficient A must vanish in a cylindrically symmetric collision, and for 
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FIG. 4: The conditional probability distribution P(if>3,3\ipi,3^ r) f° r three different impact param- 
eters 6 = 0, 7.6, 10.5 fm. The functional form of the dashed curve is given by Eq. (2.40) with fit 
parameter C = 0.53 for b = 7.6 fm and C = 0.56 for b = 10.5 fm. 



small anisotropy we have 



A oc ((e 2 > 



(2.37) 



where the double brackets denotes an event averaged e 2 Similarly dP/d(i/> 3}3 — ^/r) must 
involve even cosine terms and must be 27r/3 periodic 



dP 



d(V>3,3 - ®r) 2tt 



;i + 2A 6 cos(6(^3,3-^)) + -..) • 



(2.38) 



The relatively high fourier number n = 6 explains the smallness of the observed asymmetry, 
and Aq will be ignored from now on. The form of the conditional probability distribution can 
also be established based on general considerations. Appendix [B] uses symmetry arguments, 
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FIG. 5: The figure qualitatively describes the fluctuations associated with the Glauber model as 
illustrated in Fig. [4j When the dipole asymmetry is in plane (Position A), then the tip of triangu- 
larity is aligned with dipole asymmetry. When the dipole asymmetry is out of plane (Position B), 
the tip of the triangle is anti-aligned with the dipole asymmetry. 



a fourier expansion, and the statement that the correlation is strongest when the triangle 
and dipole angles are aligned at 7r/2 out of plane, to establish a three parameter functional 
form which describes the correlations fairly well 

^3,3^1,3, = ^ [l - 2 (Bo - 2B 2 cos(2^ lj3 - 2* R ) ) cos(3^ 3 ,3 - <P* ~ 2* *)1 , (2.39) 
where 

0* = Vh,3 - C sin(2^ li3 - 2V R ) . (2.40) 

The signs are chosen so that Bq, B 2 , and C are positive constants in the final fits. A 
sample fit with this functional form is given in Fig. 16 of Appendix [Bj The phase angle (ff is 
illustrated by the dashed black line in Fig. |4j which is found by solving 3^3 — <p* + 2^> r = n 
for ip3,3. Although we have written the conditional probability when the reaction plane angle 
is fixed, the same arguments could have been used to determine the functional form of the 
conditional probability when the participant plane angle is fixed, i. e. 



-f > (V , 3,3|^'i,3^ r pp) = Eq. (2.39) with ty R — > \l/pp and sightly different numerical coefficients. 

(2.41) 

In the limit of small elliptic eccentricity the coefficients scale as 

B cc ((e 2 » , B 2 oc ((e 2 » 2 , C cx ((e 2 » , (2.42) 

as is shown in Appendix [B} Thus the conditional probability distribution simplifies in this 
limit to 

P{^Mi,^r) = irl 1 - 2B o cos(3^ 3 ,3 - ^1,3 - 2* fl )] , (2.43) 

which describes almost all of the essential physical features. 

The strong correlation means that if the triangular and the participant planes are known, 
then the dipole plane can be determined statistically. The probability distribution of ipi 3 
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3^3,3 — 2^pp — 7r. The points indicate the average ((cos^i^ — tp^)}) and the errorbars indicate 
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for fixed ^3 and *ffpp is approximately 

1 



2tt 



1 + 2 A cos(2^ li3 - 2^ PP ) - 2B cos(3^ 3 ,3 - ^1,3 - 2^ PP ) 



(2.44) 



Maximizing this probability we determine the most probable angle of given ^3 and 
^ pp. Neglecting the A coefficient which is significantly smaller than Bq we find 



^l?3 = 3 ^3,3 



2\I> PP - 71 . 



(2.45) 



To estimate the degree of correlation between the most probable value and ^3 we calculate 

- ((cos(^i, 3 - 3^3,3 + 2tf P p)» , (2.46) 



and illustrate the result in Fig. [6j We will use this correlation in Section VA to make a 
definite prediction for the behavior of two particle correlations with respect to the reaction 
plane. 



E. Convergence of the cumulant expansion for smooth Glauber type initial con- 
ditions 

In the previous section we introduced a cumulant expansion to characterize the response 
of the system to a set of perturbations. In this section we will study the convergence 
of the cumulant expansion. Specifically, for a smooth (optical) Glauber profile, we will 
replace the initial entropy distribution with an approximately Gaussian profile and cumulant 
corrections through forth order. The distribution of entropy in the optical Glauber model 



see Appendix A 2) is first used to calculate (r ) and (r cos 20), which determines the two 



coefficients of the Gaussian. Also the normalization (i.e. the total entropy) is the same 
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FIG. 7: (Color Online) (a) Spectra in the smooth (optical) Glauber model compared to the 
cumulant expansion. The coefficients of the Gaussian and fourth order cumulant expansions have 
been adjusted to reproduce (r 2 ), (r 2 cos 20) and (r 4 cos 20), (V 4 cos 40) respectively. The total 
entropy of the cumulant expansion is also matched to the total entropy of the glauber distribution, 
(b) Elliptic flow in the Glauber model compared to the cumulant expansion. 



between the Gaussian and the Glauber distribution. Taking the impact parameter to be 
b = 7.6 fm, Fig. [7] compares the spectra and the elliptic flow for these two distributions. In 
the next approximation, the fourth cumulants to the Gaussian are adjusted as described in 
Section II D and Appendix A 2 to reproduce the (r 4 ), (r 4 cos 20), and (r 4 cos 40) moments of 
the Glauber distribution. Fig. [7] shows that the cumulant expansion reproduces the response 
of the Glauber distribution in detail. 



III. TIME DEVELOPMENT OF THE RESPONSE 



In the previous sections we introduced a set of initial conditions with definite triangularity 
and dipole asymmetry. In this section we will show how the hydrodynamic response to these 
cumulants develops in space and time. The point here is to understand the hydrodynamics 
without the complications of freezeout and a freezeout prescription. 

To show how the dipole and triangular flow develop in time, we have generalized the 
discussion of elliptic flow given in Ref. [21]. The spatial anisotropy is characterized by the 
second moment 

(r 2 cos 20) 

e ^ — ' y 6A > 

which is a function of time in general. As the system expands, the spatial anisotropy de- 
creases and the momentum anisotropy increases. The momentum anisotropy is traditionally 
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defined with e2 P 



_ J d 2 x (T xx - Tyy) = J d 2 x (e + p)u 2 r cos 20 M 
€2p ~fd 2 x (T xx + Tw) f d 2 x [(e + + 2p] ' 1 ' } 

where u r = \J [u x ) 2 + (u y ) 2 and (j) u = tan _1 (M y /V 1 ) . This definition has its flaws since the 
numerators and denominators do not transform as components of a tensor under transverse 
boosts 1 |13j . An alternative definition is found by constructing an irreducible rank two 
tensor out of the momentum density T 0t and the flow velocity v? 

T o(i u j) _ traces = I ( T <V + T o ju i _ gijj<a u ^ . (3.3) 



Then we define 



_ / d 2 xr [T^ x u x ^ - traces] _ / d 2 x ru° [(e + p)u 2 cos 20 u ] 
62p ~ J d%r [T°°u°] ~ Jd 2 xru° [(e + p)u 2 r + e} ' 1 °' ! ' 



which is almost the same as Eq. (3.2). For the triangularity and dipole asymmetry we define 
the (reducible) third rank tensor 

r 0( VV) = ^ (T 0i u ] u l + perms) . (3.5) 

Then the traceless (or irreducible) tensor is used to define the momentum space triangular 
anisotropy 

_ / d 2 xr [T^ x u x u x ^ - traces] _ / d 2 x tu° [(e + p)ul cos 30 M ] 
C3p = J d 2 xr [T 00 u°u } J d 2 xr [T 00 u°u ] ' (3 ' 6 ^ 

and the trace is used to define momentum space dipole asymmetry 

_jd 2 xr [SflTObu'vP] Jd 2 x tu° [(e + p)v% cos 4>u] 
€lp J d 2 XT [T 00 U°U } = J d 2 XT [T 00 u°u } ' 



(3.7) 



Armed with these definitions, Fig. [8] illustrates the development of the triangular flow 
and the dipole asymmetry as a function of time. As is familiar from studies of the elliptic 
flow [13 121] , the spatial anisotropy decreases leading to a growth of the momentum space 
anisotropy. When the spatial anisotropy crosses zero, the growth of the momentum space 
anisotropy stalls. The figures also indicate that the elliptic flow, the dipole asymmetry, and 
the triangularity all develop on approximately the same time scale, r ~ (r 2 )/c s . 

Another important aspect of the flow is the transverse radial flow profile. To illustrate 
this profile we decompose the transverse flow velocity into harmonics: 

Ur(r, <p) =u° r {r) + 2u^\r) cos(0) + 2uf ) {r) cos(20) + 2uf\r) cos(30) + . . . . (3.8) 



1 This flaw is easily remedied by replacing d 2 x with the fluid three volume in the local rest frame 
d 2 xdriTu . The additional factor of u° appears naturally below. 
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FIG. 8: (Color online) (a) The spatial anisotropy of the entropy distribution e± x , €2 X , and e% x (see 
text) as a function of time for b = 7.6 fm. (b) The momentum anisotropy ei p , 62 P , and e% p (see text) 
as a function of time. The timescale in these figures should be compared to \J {r 2 )/c s ~ 5.4 fm. 

For a radially symmetric Gaussian distribution only the zero-th harmonic is present, and 
shows a linearly rising flow profile. When the elliptic deformation is added the second 
harmonic also shows a linearly rising profile. Close to the origin this behavior can be 
understood with a linearized analysis of the acoustic waves. The flow velocity in an acoustic 
analysis is the gradient of a scalar function $ which can be expanded in harmonics: 

$(r,0) = $ (o) (r) + 2$ (2) (r)cos20 + ... . (3.9) 

If $(r, (p) is an analytic function of x and y, then $( 2 ) must be quadratic for small r. 
Consequently the gradient of this function, u^\r), rises linearly at small r. Similarly, the 
triangular deformation $( 3 )(r) should be cubic at small r and the flow profile uf^ should 
be quadratic. These features are borne out by our numerical work as exhibited in Fig. |9j 
Fig. [9] also shows the flow profile of the first harmonic which results from an initial dipole 
asymmetry. The first harmonic shows a negative slope at small r followed by a quadratically 
rising profile at larger r. 

As seen from Fig. [9j the triangular and dipolar flows are biased towards the edge of the 
nucleus. In the next section we will see that due to this bias v\ and v 3 are more sensitive to 
the freezeout prescription than t>2. 



IV. PARTICLE SPECTRA: vi{pr) AND v 3 {p T ) 

Having illustrated the essential features of the hydrodynamic response, we will compute 
the particle spectra associated with these flows. As discussed above, the analysis is limited 
to a classical massless ideal gas. We will follow the time honored, but poorly motivated 
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FIG. 9: (a) The zeroth harmonic of the flow profile (see Eq. 3.8 ) for the radially symmetric Gaussian 
adopted in this work. The root mean square radius of the Gaussian is adjusted to reproduce an 
impact parameter of 7.6 fm. (b) The second harmonic of the flow profile for an elliptic perturbation, 
(c) The third harmonic of the flow profile for a triangular perturbation (d) The first harmonic of 
the flow profile for a distribution with a net dipole asymmetry. The deformations ei , ei and €3 are 
all set to 0.1. 



prescription of specifying a freezeout temperature or a freezeout entropy density. Freezeout 
temperatures in full hydrodynamic simulations with a Hadronic Resonance Gas (HRG) range 
from T = 160 MeV to T = 120 MeV [T3l ITS] . The total initial entropy and initial volume 
used in our massless ideal gas simulations were taken to be the similar to the total entropy 
and initial volume used in these full hydrodynamic simulations. The final freezeout volume 
of our massless-gas simulation is also taken to be similar to the final freezeout volume of 
these full simulations. Since entropy is conserved, this can be accomplished by adjusting the 
freezeout entropy density of the massless gas so that the entropy density equals the HRG 
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Hadron Gas Tf Q 


Hadron Gas Sf D 


Massless Gas Tf Q 


130 MeV 


4.34 fm" 3 


71 MeV 


150 MeV 


1.87 fm" 3 


96 MeV 


170 MeV 


0.77 fm" 3 


127 MeV 



TABLE I: Table of freezeout temperatures used in this work. The first two columns show freezeout 
temperatures and the corresponding entropy densities of a Hadron Resonance Gas (HRG) EOS. 
The last column shows the freezeout temperatures where the massless gas EOS used in this work 
attains the corresponding HRG entropy density. 



entropy density for a specified HRG freezeout temperature. Experience has shown that 
this is a fair way to compare different equations of state. Rather than quoting the actual 
freezeout temperature of the massless gas EOS, we will simply quote the corresponding HRG 
freezeout temperature. Thus T 170 MeV means that the actual freezeout temperature is 
such that the entropy density of a massless gas is equal to the entropy density of the HRG 
at T = 170 MeV. Table |T] shows a set of temperatures and entropy densities in a HRG model 
and the corresponding freezeout temperatures for the massless gas equation of state. 



Fig. 10 shows the momentum anisotropies as a function of time, and marks when the 
average entropy density of the system reaches a specified freezeout entropy density. Specifi- 



cally, the lines indicate when (s) in the notation of Eq. (2.6) falls below the freezeout entropy 
density indicated in Table [TJ We see that for Tf Q <=>■ 170 MeV the triangular and dipole flows 
are still developing, while for Tf Q 130 MeV the flows are almost fully developed. 
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FIG. 10: Evolution of the momentum anisotropy as a function of time at an impact parameter 
of b = 7.6 fm. The lines indicate when the average entropy density (s) falls below the freezeout 
entropy density specified by the the temperatures T 130, 150, 170 MeV. 
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Once the freezeout surface is specified the particle spectra are computed using the Cooper- 
Frye formula 

.1 \ ' r 

(4.1) 



(27r) 3 E^ = J^dV,U-P-U(X)) 



where f (E) = g exp(—E/T(X)) is the distribution function of a classical massless gas. (The 
notation here follows the review article [13] ■) Using this formula we compute the particle 
spectra and determine the associated harmonics v i, Vi and v%. For each impact parameter we 
determine the root-mean square radius and the total entropy from an optical Glauber model; 
then the Gaussian parameters are adjusted to reproduce these Glauber quantities; finally 
the simulation is run to the freezeout entropy density and the harmonics are computed. 



Fig. 11 shows how the harmonics depend on centrality and the freezeout temperature. 



Examining Fig. [IT] we see that v±, t>2 and are roughly independent of centrality. How- 
ever, it must be borne in mind that in a more complete simulation, the total entropy per 
participant is also a function of centrality and this could change the result. Here the entropy 
per participant is constant. Generally the freezeout criterion is also a function of centrality 
and this could give a substantial shape to these curves in a final simulation. Finally, when 
viscosity is included the triangularity is also a more complicated function of centrality [14J. 
This will be explored elsewhere [2T] . 

Fig. 12 shows how these harmonics depend on p?. v^ipr) an d v${pt) show a characteristic 
linear rise with p? that is a consequence of a strong radial flow [25H28] . Indeed examining 
the thermal distribution with constant temperature, we have 



e P-U/T =e -E p U T /T e p T /T u r (r,cf>) co S (4> p -<f> u ) 
~ e -E v /T e p T /Tu^ (r) cos(^-0) 



(4.2) 



X 



1 + ~^uf\r) cos(20) cos(0 p - 0) + 



2 **-,(3), 



T 



-ul 



r) cos(30) cos( 



+ 



(4.3) 



Here E p is the energy, <ft p is the particles azimuthal angle; we have adopted a non-relativistic 
approximation U T ~ 1 and assumed that the flow is approximately radial, <ft u ~ <fi. Further, 
we have neglected Ur in this discussion. Unless the momentum angle equals the spatial 
angle <p p ~ <ft, the thermal distribution is strongly suppressed by the leading Boltzmann 
factor. Thus, we arrive at a form which illustrates the linear rise of rise of v n (jpr) with p? 



e P-U/T r^ e -E v /T e p T /Tu { °\r) 



1 + 



^^(r) cos(20 p ) + ^ 3 )(r) cos(30 p ) + . . . 



(4.4) 



Examining Fig. 12, we see that t>i(pr) also displays a similar linearly rising trend at higher 
Pt after an initial dip. 
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FIG. 11: vi, V2 and V3 per unit anisotropy as a function of N part for different freezeout temperatures. 
The anisotropy parameters are all 0.1 in the actual simulations. 

V. FURTHER PREDICTIONS AND COMPARISON WITH OTHER WORKS. 
A. Further predictions 



Fig. 11 and Fig. 12 show the response of the hydrodynamic system to the deformations. 
Certainly it is premature to compare the current calculation to data. For instance, the effect 
of viscosity, resonance decays, and a lattice-based equation of state have not been included. 
These reality factors will reduce the response. Nevertheless, in order to keep the final goal 
clearly in sight, we will provisionally compare the current calculation to the Alver Roland 
fit [1] of STAR inclusive two particle correlations [29J. Further, we will suggest a number of 
additional observables which can confirm the geometric nature of the measured two particle 
correlations. 
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FIG. 12: v u (pt) for two different freezeout temperatures as described in Table [TJ The root mean 
square radius of the initial Gaussian corresponds to a radius of b = 7.6 fm 

The average over glauber configurations at fixed N pait is denoted with double brackets 
((. . .)). Then the two particle angular correlation function can be expanded in a Fourier 
series: 



dK 



pairs, a/3 



— n- 



(5.1) 



The particle labels a and (3 could denote distinct particle types or px bins for example. 
Following Alver and Roland we will approximate the two particle correlation with the dis- 
connected component. The yield of particle type a for a fixed Glauber configuration is 



dN a N a 



2tt 



d(p c 



1 + 2— — ei cos( 
ei 



^i, 3 ) + 2^e 2 cos(20 a - 2^ PP ) 



+ 2^e 3 cos(30 a - 3^3,3) , (5.2) 
£3 



where we have assumed that the response is linearly proportional to the deformation. Then 
the two particle correlation function is approximated as 



diVp a i rsct| g 



dN dN 



NgN/3 

(2vr) 2 



1 + E 2 



^naVn/3 



{& cos(n(0 a - <j> p )) 



(5.3) 



Here and below we have tacitly assumed that the multiplicity fluctuations at fixed iVp art are 

N a Np ((N a N p )) ((el)) 



negligible. If this is not the case then one has the following replacements in Eq. (5.3) 



{{N a Np)) • 



(5.4) 
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FIG. 13: Fourier components of the two particle correlation function as a function of iVp a rt relative 
to the quadrapole component, (a) The triangularity component compared to the Alver Roland fit 
[1] of STAR inclusive two particle correlation functions |29j . (b) The dipole component relative to 
the quadrapole component; note that the scales differ between(a) and(b). 



Given the parameterizations in Eq. (5.1) and Eq. (5.3), the response functions in Fig. 11 



make a definite prediction for the different Fourier components V n A- The elliptic flow is 
too large in the ideal massless gas model considered here. We will therefore simply plot 
the ratios of the different fourier components as was done in the Alver and Roland paper. 
Using the response functions in Fig. 11 and the Glauber estimates for ((e§)) / ((e|)), Fig. |l3[a) 



compares the strength of the triangular component to the quadrapole component. The ideal 
hydrodynamic prediction (with a massless ideal gas EOS) is generally too large and fairly 
sensitive to the freezeout temperature. This sensitivity reflects the fact that the triangular 
flow develops further towards the edge of the nucleus. Fig. [l3](b) compares the dipole 
component to the quadrapole component. The dipole component is a factor of eight smaller 
than the quadrapole component. This is a reflection of the fact that e\ is small, and the 
fact that vi(pr)/€i is positive and negative. The dipolar flow is also sensitive to the details 
of freezeout. 

Next we wish to determine the general form of the two particle correlation function with 
respect to the participant plane ^pp 



dA^p a i rS)Q( g 
d0id02 



dN n 



dN, 



d(0 c 



pp) 



pp) 



(5.5) 



Inserting Eq. (5.2) into Eq. (5.5) and averaging over glauber configurations several sev- 



eral terms appear. In Section [B] we identified the principle correlations that exist be- 
tween the angles ^13, "03,3 an d ^pp- Namely, the only significant fourier expectation 



values are (cos(2?/> lj3 — 2^/ PP )) (as determined by the coefficient A in Eq. (2.36)), and 



(cos(^i j3 — 3^ 3j 3 + 2\I/pp)) (as determined by the coefficient B in Eq. (2.43)). With the 
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assumption that these are the only significant fourier expectation values at third order, the 
form of the two particle correlation function with respect to participant plane becomes: 



diVp airSja/3 \\ N a Np r ^ / v na v nP 



d0 Q d0 /3 // (2vr)2 L 4^ V e 



el)) cos(nd) n — n 



2 / « n 
+ 2^((e 2 » cos(20 Q - 2^ PP ) 
+ 2 V ^^((e 2 2 )) cos(20 Q + 2^ - 4* PP ) 

e 2 

+ 2^^((e?cos(2^i,3-2*pp))) cos(0 a + 0^ - 2^ PP ) 
e i 

+ 2^^(( eie3 cos(^i,3 - 3^3,3 + 2V PP ))) cos(0 a - 30^ + 2^ PP ) 

ClC3 



(5.6) 



The symmetrization with respect to a and /3 applies to all terms in this expression which 
are not already symmetric, e.g. cos(20 a — 2\I/ pp). We will discuss this expression line by 
line. The first three lines are not particularly novel: The first line is independent of the 
reaction plane angle ^ pp . The next two lines reflect the underlying elliptic flow and would 
normally be subtracted in a flow subtracted two particle correlation function. 

The fourth line contains the first novel feature. This term arises because the dipole 
asymmetry is preferentially oriented out plane, leading to a net v% out of plane. Fig. |i~4|(a) 
shows the correlation function (cos(0 a + 0^ — 2ty PP )) as a function of centrality. Recently, 
the STAR collaboration measured a similar expectation value, but divided correlation func- 
tion into the different possible charge components (i.e.++, +- , — ) in order to investigate 
the possibility of local parity violation in heavy ion collisions [30H32] . Fig. [I4|(b) shows the 
measured STAR correlations. The measured correlation is the same order of magnitude as 
the out of plane flow found in this work. However many aspects of the out of plane dipole 
flow (e.g. the pt dependence and most importantly the charge dependence) do not agree 
with the measured correlation. Thus the STAR measurements can constraint the geometric 
fluctuations reported here. This will be investigated in future work. 

A second novel feature expressed by the two particle correlation function with respect 



to reaction plane is recorded by the 5th line of Eq. (5.5). It shows that hydrodynamics, 
together with the geometric fluctuations of the Glauber model makes a definite prediction 
for the angular correlation 

((cos(0 a -30^ + 2^ PP ))). (5.7) 

Taking a to label all the particles in a definite pt bin and /3 all the particles, this definite 
prediction reads 

((cos(0 a - 30/3 + 2tfpp)» = ^^-((eieacos^ - 3^ 3 ,3 + 2tf PP )» . (5.8) 

ei £3 



This result is illustrated in Fig. 15 and is based on the Glauber analysis in Fig. p\ and the 



response functions calculated in Fig. 12 Another way to probe this same correlation is the 



following. Experimentally, the participant plane \I/pp is traditionally estimated by using 
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FIG. 14: (a) The expectation value (cos(^> Q + <pp — as predicted by hydrodynamics, 

where a and j3 label all particles, (b) The charge asymmetry with respect to reaction plane 
(cos((p a + (ftp — 2\Pfl)) as measured by the STAR collaboration [30l[31] . Here a and f3 label ++, H — , 
or . The hydrodynamic prediction does not explain the charge asymmetry. 



the standard Q vector method, or the Yang-Lee zero generalization of this idea [33]. These 
same methods can be used to determine the triangularity event plane ^3,3 without significant 
modifications [31]. The strong correlation between the dipole, triangular, and participant 
planes implies that the knowledge of ip3 3 an d ^pp determines the dipole event plane ipi,3 



at least statistically. The most probable orientation is given by Eq. (2.45) and is repeated 
here for convenience 

Thus, the v\ associated with the dipole asymmetry can be determined by measuring the 
expectation value 

<(cos(0-^3 P )>>- (5.9) 

Essentially this correlation is a v± with an extra twist to take out the shifting orientations 
of the dipole and triangular event planes - see Fig. [5] 



B. Discussion and comparison with other works 

We hope that the cumulant expansion presented in Section [TT] organizes and formalizes 
the study of fluctuations in heavy ion collisions. The convergence of the cumulant expansion 
is really quite good as illustrated in Fig. [7] At third order in the cumulant expansion there 
are two additional terms, the triangularity (r 3 cos 3(0 — ^33)), an d the dipole asymmetry 
(r 3 cos(0 - ^1,3)). 
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FIG. 15: A hydro dynamic prediction for the expectation value (cos(0 a — 30^ + 2^pp)) which 
reflects the correlation between the dipole, triangular, and elliptic event planes. Here a labels all 
particles in a given pr bin and (5 labels all particles. 

Our numerical results for the triangularity ^3/63 are similar to recently reported results 
[lU [15]. However, V3 (and v±) is significantly more sensitive to the freezeout dynamics. 
To understand this we studied the space time development of the triangularity (and dipole 
asymmetry) in Figs. [8] and [9] These figures indicate that the triangular flow develops on 
the same time scale as the elliptic flow. (A similar conclusion for the triangular flow was 
reached in Fig. 3 of Ref. [2] based on kinetic theory calculations.) However, there is an 
important difference between the elliptic flow and the dipole and triangular flows which has 
not been fully clarified previously. Specifically, the dipole and triangular moments of the 
transverse flow grow quadratically with radius, oc r 2 , rather than linearly as is the case 
with elliptic flow, u T ! oc r. Consequently, edge effects can significantly reduce the dipole 
and triangular flows. Increasing the freezeout temperature cuts on the exterior region of 
the flow profile, and therefore v 1 and v% are more sensitive to the precise freezeout criterion 



(see Figs. 11 and 12). This unfortunate result may limit the usefulness of the dipole and 
triangular flows in determining the shear viscosity of the quark gluon plasma. Indeed the 
strong reduction of the v 3 due to the shear viscosity [HI EES] is presumably largely due to 
the shear viscosity below T c , though this conclusion requires further investigation. 

We also investigated the dipole asymmetry, (r 3 cos(0 — ^1,3))- The dipole asymmetry 
appears to the same order in the gradient expansion and has not been studied previously to 
our knowledge. The dipole asymmetry is generally smaller than the triangularity since e\ is 
comparatively small. However, V\je\ is only marginally smaller than ^2/^2 and v^/e^. In non- 
central collisions the dipole asymmetry is strongly correlated with the triangularity and the 
reaction plane as is illustrated in Fig. [4] and explained in Fig. [5j We find that in non-central 
collisions the dipole asymmetry is preferentially out of plane leading to a Vi out of plane. 
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The size of the observed correlation is somewhat smaller than the observed correlations 
measured by the STAR collaboration and does not explain the charge asymmetry. 

Finally, we noted that the strong correlation between the dipole asymmetry and the trian- 
gularity can be measured experimentally by measuring two particle correlations with respect 
to reaction plane. The final result is a hydrodynamic prediction for a curious correlator 

((cos(0 a -30^ + 2^ PP ))), (5.10) 

which is shown in Fig. [T5j This average is similar to averages used to investigate the 
Chiral Magnetic Effect (CME) and is no more difficult to measure. The hydrodynamic 



prediction for Eq. (5.10) is several times larger than the correlation currently measured 



by the STAR collaboration, (cos(0 a + — 2^pp)). Thus, the proposed measurement is 



feasible and important. If the predictions of Fig. [15] are confirmed it would validate the 
hydrodynamic and geometric nature of the measured two particle correlations. Further, 
given the off-diagonal nature of the proposed measurement, it will be difficult to reproduce 
this correlation with other mechanisms. 

The current study neglected the effects of shear viscosity and resonance decays and used 
an ideal gas rather than a lattice based equation of state. Incorporating these important 
corrections is left for future work. 
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Appendix A: Details of the cumulant expansion and initial conditions 
1. Formal expansion 

Our goal is to determine the cumulants of the underlying distribution p(x) and to de- 
compose these cumulants into irreducible tensors with respect to rotations around the z 
axis. 

First we expand p(x) and its Fourier transform p(k) in a fourier series 

p(x) = p(r, 0) =p Q (r) + 2 p c n {r) cos(n0) + 2 ^ p s n {r) sin(n0) , (Al) 

n=l n=l 

oo 

p{k) = p{k, <f) k ) =p {k) + Pnik) sin(n0 fe ) + 2 ^ p s n (k) sin(n0 fe ) , (A2) 

n=l n=l 

where r,<ft,k,<f)k are the magnitudes and azimuthal angles of x and k respectively. The 
relation between the pn S {k) and p^' s (r) is established by substituting the identity 

oo 

e ih x = J (kr) + 2^2 i n Jni kr ) c os(0 - 0ft) ( A 3) 

n=l 
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into the Fourier transform (Eq. ( 2.3[ )) and using elementary manipulations to obtain 



p c /(k) = 2tt / rdri n J n (kr)p c n s (r). 



(A4) 

Similarly, the generating function of cumulants is also given by a fourier series 

W{k) = W (k) + 2 W n(k) cos(n0 fc ) + 2 ^ ^(fc) sin(n0 fc ) , (A5) 



and each W% s {k) is expanded in k as described by equations Eqs. 2.16 and 2.22 Then we 
can expand both sides of the defining relation 



exp(W(k)) = p{k) 



(A6) 



in a series expressions of the form k m cos(n0k) and k m sin(n0fc). In developing this expansion 
we use the series expansion of the Bessel function 



oo / [ 2\k 



k=0 



k\T(u + k + 1) 



(A7) 



and the series expansion of W^ ,s (k). Comparing idential powers of k m cos(n<f)k) and 
k m sin(n0fc) we determine the W% s m in terms of the moments of the underlying distribution. 
Through fifth order inclusive this comparison yields the following relations: 



0-th harmonic: 



2nd harmonic: 



o.i 



( r 4 )-2(r 2 ) 2 -(r 2 cos20)' 



(A8) 
(A9) 
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r 2 cos : 



>>]. 



[<r 4 cos 20) - 3 (r 2 ) <r 2 cos 20)] 



\{r sin 2 



(A10) 
(All) 
(A12) 



4th harmonic: 



(r 4 cos 40) -3<r 2 cos(20)) 

W, 



; ^-[<r 4 sin40>] 



(A13) 
(A14) 
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1st harmonic: 



8 

W7,s=jj [<r 3 sin(0)>] , 
5 r 

(r 5 cos(0)) — 6 (r 2 ) (r^ cos0) 



,5 



3rd harmonic: 
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16 



- yir 2 cos 20) <r 3 cos 30) + 3 (r 3 cos 0) (r 2 cos 20) 
(r 3 sin(0)) — 6 (r 2 ) (r 3 sin0) 



r 2 cos 20\ (r 3 sin 36) — 3 (r 3 sin d>) (r 2 cos 20*1 



3,3 
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[<r 3 cos(30)>] , 
[<r 3 sin(30)>] , 



W. 



3,5 



^3,5 



32 L 
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32 



(r 5 cos 30) — 4 (r 2 ) (r 3 cos 30) — 6 (r 3 cos 0) (r 2 cos ! 



(r 5 sin 30) - 4 (r 2 ) (r 3 sin 30) - 6 (r 3 sin 0) (r 2 cos 2 



5th harmonic: 



i 

32 
1 

32 



(r 5 cos(50)) - 10 (r 2 cos 20) (r 3 cosd 



(r 5 sin(50)) - 10 (r 2 cos 20) (r 3 sin 3. 



Each coefficient has a simple interpretation. For instance, Wq^ 



(A15) 
(A16) 

(A17) 

(A18) 

(A19) 
(A20) 
(A21) 
(A22) 

(A23) 
(A24) 

is simply the 



root mean square radius of the Gaussian. To classify corrections to the Gaussian, one should 



examine the difference between (r 4 ) and (r 2 ) 2 ; W 0j4: ~ | (r 4 ) — 2 (r 2 y is the required 

difference. The underlined terms (i.e. (r 2 cos20) 2 in the case W 0:4 ) are of suppressed by a 
power of e 2 and are therefore generally unimportant except in very peripheral collisions. 
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2. Fourier transform and regulating the cumulant expansion 



After specifying the cumulants, the distribution is Fourier transformed to determine the 
initial entropy density in coordiante space. For simplicity we will discuss only a spherically 
symmetric Gaussian deformed by a small definite triangularity, W£ 3 . In this case the Fourier 
transform of a distribution with Wo j2 and small W£ 3 , 



p(x) 



d 2 k 

(27T) 2 



-ikx 



e 2 



W ,2 



l + ^3 C , 3 W 3 cos30 fc + 



(A25) 
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yields with the definition of €3 in Eq. (2.30) 



p(x) 



1 + 



(r 3 ) e 3 
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r 2 
<r2> 



7r (r 2 ) 



(A26) 



where (r 3 ) = 3^/4 (r 2 ) 3/2 . At large enough radius the correction term becomes large 
compared to the leading Gaussian. To regulate this term we replace the whole correction 
(= X) with 

X -> Ctanh(X/C), (A27) 

where C = 0.95. We have checked that the results are independent of the precise value of 
the constant C. The regulator here is not perfect as it (weakly) mixes different terms in the 
fourier expansion, but we have found this to be unimportant from a practical perspective, 
i.e. the v% produced by this regulated €3 distribution is small. Anot her complication is 
that the input parameter e 3 nput in the regulated version of Eq. (A26) does not actually 



equal the "true" 63 of the initial distribution. In all figures we have divided by the true 
e 3 = — (r 3 cos(30)) rather than the input parameter e 



input 



We can now specify precisely the initial conditions that are used for Fig. [TT] and other 
results. At a given impact parameter we use the optical glauber model to calculate the 
distribution of participants the transverse plane with <jnn = 40 mb. In a traditional hydro- 
dynamic simulation (such labeled by the "Glauber" curves in Fig. [7]) the entropy density at 
an initial time r n = 1 fm is 



C„ dN n 



r dxdy' 



(A28) 



where — - 



is the number of participants per unit area. The value C s = 15.9 closely cor- 
responds to the results of full hydrodynamic simulations [251 133 ES] The equation of state 
that is used in this work is a classical massless ideal gas V = 1/3 e. The relation between 
the temperature and energy density is e/T 4 ~ 12.2 which is the value for a two flavor ideal 
quark-gluon plasma. In the current simulations we calculate the total entropy and average 
squared radius (r 2 ) for glauber distribution. We then take a deformation e 3 ~ 0.1, and use 
these parameters to initialize the regulated Gaussian described by Eq. ( A26[ ) and Eq. (A27). 
Finally, the simulation is run and the spectra are calculated leading to Fig. [2] 



Appendix B: Correlations in the Glauber model 



The goal of this appendix is to motivate Eq. (2.39). A given distribution of participants 
is first characterized by the participant plane ^fpp = ip2,2 and we will assume that €2 is 
small. Then the probability distribution for ^13 for fixed ^>pp is given by Eq. (2.36). For 
fixed typp and ipi,3 the probability distribution for ^3 must be 27r/3 periodic. Measuring 
all angles with respect to participant plane and keeping only the first non-trivial term in the 
Fourier series we have 



1 

27 



1 + 2B cos (3(^3,3 - *pp) - (</>* - *pp)) 



(Bl) 



28 




%I2 n 3n/2 2n 

V1, 3 -^R 



FIG. 16: A fit based on Eq. |2.39| to the the Glauber data exhibited in Fig. [4} The parameters 
are B = 0.277(2), B 2 = 0.029(1), and C = 0.532(7). The normalization (i.e. the color scale) is 
arbitrary, but is the same as in Fig. [4j 



The amplitude B and phase 0* are functions of ^3 — ^pp. 

The amplitude B and the phase derivative can be expanded in a Fourier series 

B =B + 2B 2 cos (2^1,3 - 2^ PP ) , 
:Co + 2C 2 cos(2^ li 3-2*p P ) . 



dtp 



1,3 



(B2) 
(B3) 



As the ipi j3 increases by 2n, the phase <fi* must change by a multiple of 2tt to leave the 
conditional probability distribution invariant. The simplest possibility which qualitatively 
describes the trends illustrated in Fig. [|and Fig. [5] is to take Co = 1. In a general fourier 
series of two variables other possibilities would be allowed, e.g. Cq = 3. However such 



correlations turn out to be small in the Glauber model. Integrating Eq. (B3) we find 



+ C 2 sin(2'0i ) 3 - 2^ PP ) + const 



(B4) 



The constant required to reproduce Fig. [5] is it. The combination of Eqs. Bl, B2, and B4 



leads to the parameterization quoted in Eq. (2.39). In Eq. (2.39) we absorbed the constant 



phase 7r into the leading minus sign of B and B 2 and changed the sign of C 2 so that all 



coefficients are positive in the final fit. Fig. 16 shows a fit to the Monte Carlo Glauber 
shown in Fig. [5] at b = 7.6 fm using this parameterization. The fit does capture most of the 
essential features, but fails to reproduce the sharpness of the correlation band. 

Finally, we can estimate the scaling of these coefficients with the average elliptic eccen- 
tricity ((e 2 )). In a central collision B(ipi 3 , ^pp) must vanish. This can be understood by 
examining Fig. [5] and recognizing that in a central collision there is no distinguishable dif- 
ference between Position A and Position B. The coefficient of cos(3'03 j 3 — <p* — ^pp) (i.e. 
B) describes how phase between the triangular and the dipole planes changes from Position 
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A to Position B. This coefficient must vanish in central collisions where Position A and 
Position B are identical. Finally the coefficients B2 and C2 reflect the almond shape and 
must involve an additional power of ((62)) relative to Co and Bq. With these remarks we 
arrive at the scalings given in Eq. (2.42). 
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